Temporal Cox model
Hening Cui
5/9/2023
library(Seurat)
library(spacexr)
library(tidyverse)
library(spatstat)
library(nlme)
library(lgcp)
library(CAinterprTools)
library(sf)
library(concaveman)
library(ggplot2)loading dataset
load(file= "/Users/cuihening/Desktop/shikun/seurat/robjs/chicken_visium.4.Robj")TM <- FetchData(object=chicken_visium,vars=c('TMSB4X','ident', 'orig.ident'))
tmsb4x=cbind(cell = rownames(TM), TM) %>%
dplyr::rename(time = orig.ident)
tmsb4x$is2=as.factor(ifelse(tmsb4x$ident==2,1,0))The function to get coordinate
getcoor=function(slice, time){
coords <- GetTissueCoordinates(chicken_visium, image= slice)
colnames(coords) <- c("x", "y")
coords <- cbind(cell = rownames(coords),coords)
coords$time = time
return(coords)
}coords_d4 = getcoor('slice1','D4')
coords_d10 = getcoor('slice1_D10-C1','D10')
coords_d7 = getcoor('slice1_D7-B1','D7')
coords_d14 = getcoor('slice1_D14-D1','D14')Function for constuct ppp
buildppp=function(dataset,markers){
X <- dataset$x
Y <- dataset$y
xrange <- range(X, na.rm=T)
print(xrange)
yrange <- range(Y, na.rm=T)
print(yrange)
ppppro = ppp(X,Y,xrange,yrange,marks=markers)
return(ppppro)
}EDA
D4
D4 = tmsb4x %>%
filter(time=='D4') %>%
merge(coords_d4, by="cell") %>%
filter(x >= 385 & y>=267 ) %>%
rename(Y=x, X=y) %>%
rename(x=X, y=Y) %>%
mutate(x=x-138, y=y-208) %>%
mutate(time=4)
summary(D4)## cell TMSB4X ident time.x is2
## Length:141 Min. :0.04300 6 :66 Length:141 0:119
## Class :character 1st Qu.:0.05200 4 :27 Class :character 1: 22
## Mode :character Median :0.05554 2 :22 Mode :character
## Mean :0.05688 8 :21
## 3rd Qu.:0.05953 5 : 5
## Max. :0.07820 0 : 0
## (Other): 0
## y x time.y time
## Min. :182.6 Min. :142.5 Length:141 Min. :4
## 1st Qu.:202.2 1st Qu.:180.0 Class :character 1st Qu.:4
## Median :221.8 Median :206.3 Mode :character Median :4
## Mean :222.8 Mean :201.6 Mean :4
## 3rd Qu.:241.3 3rd Qu.:225.1 3rd Qu.:4
## Max. :274.0 Max. :247.5 Max. :4
## xym= D4 %>% select(c("x","y"))
pts <- st_as_sf(xym, coords=c('x','y'))
conc <- concaveman(pts)
ggplot() +
geom_sf(data = conc, fill = NA)
d4_owin= as.owin(as_Spatial(conc))
plot(d4_owin)
D4_samppp = ppp(D4$x,D4$y, c(0, 450),c(0,550), marks=D4$is2)D4_samppp = ppp(D4$x,D4$y, window = d4_owin, marks=D4$is2)## Warning: 19 points were rejected as lying outside the specified windowD4_ty2=D4 %>%
filter(is2==1)
D4_ty2p=ppp(D4_ty2$x,D4_ty2$y, c(0, 450),c(0,550))D4_ty2_owin=ppp(D4_ty2$x,D4_ty2$y, window = d4_owin)## Warning: 4 points were rejected as lying outside the specified windowpar(mfrow=c(1,2))
plot(D4_samppp, pch = 20,show.window=FALSE,cols=c("blue", "green"), legend=FALSE, main="D4 point",markscale = 0.1)## Warning in plot.ppp(D4_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 19 illegal points also plottedplot(density(D4_ty2_owin), main="D4 density")
pixel=pixellate.ppp(D4_ty2_owin,fractional=TRUE)
plot(pixel)
plot(as.im(density(D4_ty2_owin), fractional=TRUE))
par(mfrow=c(1,2))
plot(D4_samppp, pch = 20,show.window=FALSE,cols=c("blue", "green"), legend=FALSE, main="D4 point",markscale = 0.1)## Warning in plot.ppp(D4_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 19 illegal points also plottedplot(density(D4_ty2p), main="D4 density")
D4_fit = kppm(D4_ty2p ~1, "LGCP")
D4_fit## Stationary Cox point process model
## Fitted to point pattern dataset 'D4_ty2p'
## Fitted by minimum contrast
## Summary statistic: K-function
##
## Uniform intensity: 8.888889e-05
##
## Cox model: log-Gaussian Cox process
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 5.041248 62.154818
## Fitted mean of log of random intensity: -11.84875D4_fit = kppm(D4_ty2_owin ~1, "LGCP")
D4_fit## Stationary Cox point process model
## Fitted to point pattern dataset 'D4_ty2_owin'
## Fitted by minimum contrast
## Summary statistic: K-function
##
## Uniform intensity: 0.002988242
##
## Cox model: log-Gaussian Cox process
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 1.036636e+00 3.273900e+07
## Fitted mean of log of random intensity: -6.331388plot(predict.dppm(D4_fit))
D7
D7 = tmsb4x %>%
filter(time=='D7') %>%
merge(coords_d7, by="cell") %>%
filter(x >= 255 & y>=267 ) %>%
rename(Y=x, X=y) %>%
rename(x=X, y=Y) %>%
mutate(x=x-150, y=y-152) %>%
mutate(time=7)
summary(D7)## cell TMSB4X ident time.x is2
## Length:505 Min. :0.00000 1 :166 Length:505 0:438
## Class :character 1st Qu.:0.05050 5 : 83 Class :character 1: 67
## Mode :character Median :0.05403 2 : 67 Mode :character
## Mean :0.05413 4 : 53
## 3rd Qu.:0.05735 8 : 48
## Max. :0.07924 7 : 35
## (Other): 53
## y x time.y time
## Min. :107.8 Min. :131.1 Length:505 Min. :7
## 1st Qu.:179.5 1st Qu.:172.2 Class :character 1st Qu.:7
## Median :225.1 Median :198.6 Mode :character Median :7
## Mean :220.9 Mean :200.3 Mean :7
## 3rd Qu.:264.3 3rd Qu.:228.6 3rd Qu.:7
## Max. :323.1 Max. :281.1 Max. :7
## xym= D7 %>% select(c("x","y"))
pts <- st_as_sf(xym, coords=c('x','y'))
conc <- concaveman(pts)
ggplot() +
geom_sf(data = conc, fill = NA)
d7_owin= as.owin(as_Spatial(conc))
plot(d7_owin)
D7_samppp = ppp(D7$x,D7$y, window = d7_owin, marks=D7$is2)## Warning: 34 points were rejected as lying outside the specified windowD7_ty2=D7 %>%
filter(is2==1)
D7_ty2p=ppp(D7_ty2$x,D7_ty2$y, window = d7_owin)## Warning: 3 points were rejected as lying outside the specified windowpar(mfrow=c(1,2))
plot(D7_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D7 point")## Warning in plot.ppp(D7_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 34 illegal points also plottedplot(density(D7_ty2p), main="D7 density")
D7_fit = kppm(D7_ty2p ~1, "LGCP")
D7_fit## Stationary Cox point process model
## Fitted to point pattern dataset 'D7_ty2p'
## Fitted by minimum contrast
## Summary statistic: K-function
##
## Uniform intensity: 0.002783736
##
## Cox model: log-Gaussian Cox process
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 8.909741e-01 1.332978e+08
## Fitted mean of log of random intensity: -6.329448D10
D10 = tmsb4x %>%
filter(time=='D10') %>%
merge(coords_d10, by="cell") %>%
filter(x >= 193 & y>=320) %>%
rename(Y=x, X=y) %>%
rename(x=X, y=Y) %>%
mutate(x=x-217, y=y-128) %>%
mutate(time=10)
summary(D10)## cell TMSB4X ident time.x is2
## Length:888 Min. :0.03259 3 :163 Length:888 0:760
## Class :character 1st Qu.:0.05015 2 :128 Class :character 1:128
## Mode :character Median :0.05414 1 :123 Mode :character
## Mean :0.05418 4 :118
## 3rd Qu.:0.05806 0 :109
## Max. :0.07277 7 :108
## (Other):139
## y x time.y time
## Min. : 72.12 Min. :103.4 Length:888 Min. :10
## 1st Qu.:169.94 1st Qu.:156.0 Class :character 1st Qu.:10
## Median :215.79 Median :197.4 Mode :character Median :10
## Mean :220.37 Mean :200.8 Mean :10
## 3rd Qu.:267.97 3rd Qu.:238.5 3rd Qu.:10
## Max. :378.87 Max. :321.0 Max. :10
## xym= D10 %>% select(c("x","y"))
pts <- st_as_sf(xym, coords=c('x','y'))
conc <- concaveman(pts)
ggplot() +
geom_sf(data = conc, fill = NA)
d10_owin= as.owin(as_Spatial(conc))
plot(d10_owin)
D10_samppp =ppp(D10$x,D10$y, window = d10_owin, marks=D10$is2)## Warning: 46 points were rejected as lying outside the specified windowD10_ty2=D10 %>%
filter(is2==1)
D10_ty2p=ppp(D10_ty2$x,D10_ty2$y, window = d10_owin)## Warning: 2 points were rejected as lying outside the specified windowpar(mfrow=c(1,2))
plot(D10_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D10 point")## Warning in plot.ppp(D10_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 46 illegal points also plottedplot(density(D10_ty2p), main="D10 density")
D10_fit = kppm(D10_ty2p ~1, "LGCP")
D10_fit## Stationary Cox point process model
## Fitted to point pattern dataset 'D10_ty2p'
## Fitted by minimum contrast
## Summary statistic: K-function
##
## Uniform intensity: 0.002817653
##
## Cox model: log-Gaussian Cox process
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 1.000520e+00 1.545034e+06
## Fitted mean of log of random intensity: -6.372111D14
D14 = tmsb4x %>%
filter(time=='D14') %>%
merge(coords_d14, by="cell") %>%
rename(Y=x, X=y) %>%
rename(x=X, y=Y) %>%
mutate(x=x-144, y=y-78) %>%
mutate(time=14)
summary(D14)## cell TMSB4X ident time.x is2
## Length:1967 Min. :0.03193 0 :854 Length:1967 0:1643
## Class :character 1st Qu.:0.05184 2 :324 Class :character 1: 324
## Mode :character Median :0.05530 3 :175 Mode :character
## Mean :0.05509 4 :159
## 3rd Qu.:0.05871 8 :138
## Max. :0.07352 7 :128
## (Other):189
## y x time.y time
## Min. : 9.428 Min. : 58.58 Length:1967 Min. :14
## 1st Qu.:139.719 1st Qu.:137.43 Class :character 1st Qu.:14
## Median :224.323 Median :200.63 Mode :character Median :14
## Mean :220.551 Mean :200.92 Mean :14
## 3rd Qu.:302.755 3rd Qu.:261.00 3rd Qu.:14
## Max. :400.476 Max. :361.98 Max. :14
## xym= D14 %>% select(c("x","y"))
pts <- st_as_sf(xym, coords=c('x','y'))
conc <- concaveman(pts)
ggplot() +
geom_sf(data = conc, fill = NA)
d14_owin= as.owin(as_Spatial(conc))
plot(d14_owin)
D14_samppp = ppp(D14$x,D14$y, window = d14_owin, marks=D14$is2)## Warning: 64 points were rejected as lying outside the specified windowD14_ty2=D14 %>%
filter(is2==1)
D14_ty2p=ppp(D14_ty2$x,D14_ty2$y, window = d14_owin)## Warning: 11 points were rejected as lying outside the specified windowpar(mfrow=c(1,2))
plot(D14_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D14 point")## Warning in plot.ppp(D14_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 64 illegal points also plottedplot(density(D14_ty2p), main="D14 density")
D14_fit = kppm(D14_ty2p ~1, "LGCP")
D14_fit## Stationary Cox point process model
## Fitted to point pattern dataset 'D14_ty2p'
## Fitted by minimum contrast
## Summary statistic: K-function
##
## Uniform intensity: 0.00337893
##
## Cox model: log-Gaussian Cox process
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 1.216153 4084.813883
## Fitted mean of log of random intensity: -6.298273Spatstat replicate model
Create the list of all point
time_ty2=solist(D4_ty2_owin,D7_ty2p, D10_ty2p, D14_ty2p)
time_ty2## List of point patterns
##
## Component 1:
## Planar point pattern: 18 points
## window: polygonal boundary
## enclosing rectangle: [142.5205, 247.5368] x [182.562, 274.0155] units
## *** 4 illegal points stored in attr(,"rejects") ***
##
## Component 2:
## Planar point pattern: 64 points
## window: polygonal boundary
## enclosing rectangle: [131.0609, 281.1149] x [107.7713, 323.0721] units
## *** 3 illegal points stored in attr(,"rejects") ***
##
## Component 3:
## Planar point pattern: 126 points
## window: polygonal boundary
## enclosing rectangle: [103.3981, 320.9683] x [72.1171, 378.8714] units
## *** 2 illegal points stored in attr(,"rejects") ***
##
## Component 4:
## Planar point pattern: 313 points
## window: polygonal boundary
## enclosing rectangle: [58.5757, 361.9798] x [9.428, 400.4762] units
## *** 11 illegal points stored in attr(,"rejects") ***sapply(time_ty2, npoints)## [1] 18 64 126 313plot.solist(time_ty2, main="The ty2 cell", main.panel=letters[1:4], legend=FALSE, equal.scales = TRUE, halign=TRUE, valign=TRUE)
Turn list into hyperframe
hyoertime <- hyperframe(cell=time_ty2)
hyoertime$time=c(4,7,10,14)plot(hyoertime, quote(plot(density(cell), main=time)))
try=mppm(cell~time, hyoertime)test the model
cdf.test.mppm(try, "y")##
## Spatial Kolmogorov-Smirnov test
##
## data: predicted cdf of covariate '"y"' evaluated at data points of 'try'
## D = 0.59049, p-value < 2.2e-16
## alternative hypothesis: two-sidedplot.mppm(try, main="The ty2 cell", cif=TRUE,
how="image")## Warning: glm.fit: algorithm did not converge
lurking(try, covariate = expression(time))## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not converge
sub_time=subfits(try)## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: algorithm did not convergefitted_intensities <- predict.mppm(try)fitted_intensities2 <- predict.mppm(try,locations = D4_ty2_owin)
fitted_intensities2## Hyperframe:
## trend cif
## 1 (im) (im)
## 2 (im) (im)
## 3 (im) (im)
## 4 (im) (im)K_functions4 <- Kinhom(D4_ty2_owin, ratio= TRUE)
K_functions7 <- Kinhom(D7_ty2p, ratio= TRUE)
K_functions10 <- Kinhom(D10_ty2p, ratio= TRUE)
K_functions14 <- Kinhom(D14_ty2p, ratio= TRUE)Kcomb=pool(K_functions4,K_functions7,K_functions10,K_functions14)try=lgcp.estK(Kcomb)try## Minimum contrast fit (object of class "minconfit")
## Model: log-Gaussian Cox process
## Covariance model: exponential
## Fitted by matching theoretical K[inhom] function to Kcomb
##
## Internal parameters fitted by minimum contrast ($par):
## sigma2 alpha
## 1.2216365337 0.0001424789
##
## Fitted covariance parameters:
## var scale
## 1.2216365337 0.0001424789
##
## Converged successfully after 57 function evaluations
##
## Starting values of parameters:
## sigma2 alpha
## 1 1
## Domain of integration: [ 0 , 22.86 ]
## Exponents: p= 2, q= 0.25Cox Spatial Temporal Model
prepare data
temp= rbind(D4, D7, D10, D14) %>%
filter(ident==2) %>%
select(x,y,time, is2)temppx=ppx(temp)
plot(temppx)
library("stpp")## Warning: package 'stpp' was built under R version 4.1.2## Loading required package: rpanel## Warning: package 'rpanel' was built under R version 4.1.2## Loading required package: tcltk## Package `rpanel', version 1.1-5: type help(rpanel) for summary information##
## Attaching package: 'rpanel'## The following object is masked from 'package:tidyr':
##
## population## Loading required package: splancs## Warning: package 'splancs' was built under R version 4.1.2## Loading required package: sp## Warning: package 'sp' was built under R version 4.1.2##
## Spatial Point Pattern Analysis Code in S-Plus
##
## Version 2 - Spatial and Space-Time analysis##
## Attaching package: 'splancs'## The following object is masked from 'package:dplyr':
##
## tribble## The following object is masked from 'package:tidyr':
##
## tribble## The following object is masked from 'package:tibble':
##
## tribbletlim=c(0,14)xyt <- stppp(list(data = temp, tlim = tlim, window = owin(c(0,450),c(0,550))))xyt <- integerise(xyt)
xyt## Space-time point pattern
## Planar point pattern: 541 points
## window: rectangle = [0, 450] x [0, 550] units
## Time Window : [ 0 , 14 ]den <- lambdaEst(xyt, axes = TRUE)
plot(den)sar <- spatialAtRisk(den)
sar## SpatialAtRisk object
## X range : [1.7578125,448.2421875]
## Y range : [2.1484375,547.8515625]
## dim : 128 x 128plot(sar)
mut <- muEst(xyt)
mut## temporalAtRisk object
## function (t)
## {
## if (length(t) > 1) {
## stop("Function only works for scalar numeric objects, ie a vector of length 1.")
## }
## if (!any(as.integer(t) == tvec)) {
## return(NA)
## }
## return(obj[which(as.integer(t) == tvec)] * scale)
## }
## <bytecode: 0x7fa201af6240>
## <environment: 0x7fa201af7468>
## attr(,"tlim")
## [1] 0 14
## attr(,"class")
## [1] "temporalAtRisk" "function"
## Time Window : [ 0 , 14 ]plot(mut)
gin <- ginhomAverage(xyt, spatial.intensity = sar, temporal.intensity = mut)##
|
| | 0%
|
|===================== | 30%
|
|========================================== | 60%
|
|======================================================================| 100%
## Returning an average of 4 curvesplot(gin)
kin <- KinhomAverage(xyt, spatial.intensity = sar, temporal.intensity = mut)##
|
| | 0%
|
|===================== | 30%
|
|========================================== | 60%
|
|======================================================================| 100%
## Returning an average of 4 curves